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Abstract 

The S0(5) D S0(3) spherical harmonics form a natural basis for expansion of nuclear collective model angular wave 
functions. They underlie the recently-proposed algebraic method for diagonalization of the nuclear collective model 
Hamiltonian in an SU(1, 1) x S0(5) basis. We present a computer code for explicit construction of the S0(5) D S0(3) 
spherical harmonics and use them to compute the Clebsch-Gordan coefficients needed for collective model calculations 
in an SO(3)-coupled basis. With these Clebsch-Gordan coefficients it becomes possible to compute the matrix elements 
of collective model observables by purely algebraic methods. 

Program Summary 

Title of program: GammaHarmonic 
Catalogue identifier: AECY_vl_0 

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AECY_vl_0 

Program available from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Operating system: Any which supports Mathematica; tested under Microsoft Windows XP and Linux 

Programming language used: Mathematica 6 

Number of bytes in distributed program, including test code and documentation: 16 037234 
Distribution format: tar.gz 

Nature of problem: Explicit construction of S0(5) D S0(3) spherical harmonics on S^. Evaluation of S0(3)-reduced 
matrix elements and S0(5) D S0(3) Clebsch-Gordan coefficients (isoscalar factors). 

Method of solution: Construction of S0(5) D S0(3) spherical harmonics by orthonormalization, obtained from a 
generating set of functions, according to the method of D. J. Rowe, P. S. Turner, and J. Repka [J. Math. Phys. 45 
(2004) 2761]. Matrix elements and Clebsch-Gordan coefficients follow by construction and integration of S0(3) scalar 
products. 

PACS: 02.20.Qs, 02.70.-c, 21.60.Ev, 21.60.Fw 
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1. INTRODUCTION 

The S0(5) D S0(3) spherical harmonics constitute the 
natural basis for the "angular" wave functions in the 
collective model of nuclear quadrupole motion [1]. The 
S0(5) D S0(3) spherical harmonics underlie the recently- 
proposed algebraic scheme [2-4] for the nuclear collective 
model. Direct products of the S0(5) D S0(3) spherical 
harmonics with appropriate optimal radial wave func- 
tions provide an SU(1, 1) x S0(5) algebraic basis [2-4] 
which allows for exceedingly efficient numerical diago- 
nalization of nuclear collective model Hamiltonians. 

For applications to transitional and deformed nuclei, 
the SU(1, 1) X S0(5) scheme reduces by orders of magni- 
tude [3] the basis size needed for convergence as compared 



to conventional diagonalization in a five-dimensional os- 
cillator [U(5) D S0(5)] basis [5-7]. Matrix elements of 
an essentially unlimited set of potential and kinetic en- 
ergy operators, often in analytic form, can easily be con- 
structed in an SU(1,1) x S0(5) basis [3, 4], in terms of 
S0(5) D S0(3) Clebsch-Gordan coefficients. The resuh- 
ing calculational scheme, the so-called algebraic collective 
model (ACM), is described in Refs. [2-4, 8]. Examples 
of physical applications may be found in Ref. [9]. 

In this article, we present a computer code for explicit 
construction of the S0(5) D S0(3) spherical harmonics 
and for using them to determine the Clebsch-Gordan co- 
efficients for coupling of symmetric irreducible represen- 
tations of S0(5) in an S0(3) basis. The construction of 
the S0(5) D S0(3) spherical harmonics is carried out by 
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orthonornialization of monomials in a set of four generat- 
ing functions, according to the method of Rowe, Turner, 
and Repka [10]. The purpose of this code is to make 
calculations using the ACM routinely possible, without 
the need to reconstruct the S0(5) machinery for each 
new application. The S0(5) D S0(3) Clebsch-Gordan 
coefficients yielded by the code are also relevant to the 
U(6) interacting boson model (IBM) [11, 12] which, in 
U(6) D U(5) D S0(5) and U(6) D S0(6) D S0(5) bases, 
is in close correspondence with the collective model in 
appropriate SU(1,1) x S0(5) bases [13]. They can, fur- 
thermore, be used as seed coefficients in the generation of 
Clebsch-Gordan coefficients for the coupling of more gen- 
eral (non-symmetric) representations of S0(5) or Sp(4) 
in an S0(3) basis [14]. 

The code accompanying this article is implemented in 
Mathematica 6 [15]. All calculations are carried out in 
exact symbolic arithmetic. A machine-readable tabula- 
tion of calculated Clebsch-Gordan coefficients is included 
along with the code in the CPC Program Library. This 
tabulation is of sufficient extent to support basic calcu- 
lations using the algebraic collective model, without the 
necessity of re-running the code. 

The basic algorithm used here is that presented in 
Ref. [10]. However, the computational techniques have 
been substantially developed. For example, integration 
via the Fourier representation of functions has been used 
to greatly increase the efficiency of the calculations and 
considerably extend the practical range of application of 
the algorithm. 

The necessary mathematical definitions for the 
S0(5) D S0(3) spherical harmonics and the general 
method for construction of the basis and evaluation of 
Clebsch-Gordan coefficients are discussed in Sec. 2. The 
more technical details of the computer implementation 
are summarized in Sec. 3. Instructions for installation 
and use of the computer code are given in Sec. 4. 



2. THE BASIC ALGORITHM 

2.1. The SO(5) D SO(3) spherical harmonics 

The S0(5) spherical harmonics are eigenfunctions of 
the Laplace-Beltrami operator A^ on the four-sphere 54, 
that is, the angular part of the Laplacian in five di- 
mensions. This operator A^ is also the second order 
Casimir invariant of S0(5). Thus, the spherical harmon- 
ics are functions of a set of coordinates on S4. Standard 
(7, il) coordinates for ^4 are reviewed below in Sec. 2.2. 
The S0(5) spherical harmonics constitute a complete or- 
thonormal basis for the space LF'{Si) of square-integrable 
functions on ^4 and transform under S0(5) rotations as 
bases for the symmetric irreps {v, 0), for u = 0, 1, . . ., of 
S0(5). 

For applications to the nuclear collective model, we 

seek spherical harmonics which have "good S0(3) angu- 
lar momentum" , that is, which also transform as bases for 



irreps, labelled by (L), of the S0(3) subalgebra of S0(5). 
The desired S0(5) D S0(3) spherical harmonics, which 
we denote by '^vaLMi'^,^), thus reduce the subalgebra 
chain 



SO(5)dSO(3)dSO(2), 

V a L M 



(1) 



with the representation labels as shown. The label v 
is conventionally termed the "boson seniority" quan- 
tum number, following Racah. It is the five- dimensional 
[S0(5)] analog of the angular momentum quantum num- 
ber. The S0(5) spherical harmonics satisfy the eigen- 
value equation 



(2) 



They are also eigenfunctions of the Casimir operators of 
the other algebras in the chain (1). However, the labels 
vLM provided by these Casimir operators are insufficient 
to fully distinguish the spherical harmonics. In particu- 
lar, multiple S0(3) representations of the same L may 
occur within a given S0(5) representation. The "missing 
label" is provided by a multiplicity index a. The branch- 
ing rule for S0(3) irreps occurring within an S0(5) ir- 
rep is well known [16, 17], and the multiplicity is given 
by (Al). 

The S0(5) D S0(3) spherical harmonics have physical 
significance as the angular, i.e., (7, fl), wave functions for 
the nuclear collective model, for the case in which which 
the collective potential is S0(5) invariant [18, 19]. Recall 
the quantum mechanics of a particle in three-dimensional 
Euclidean space, subject to a central force, i.e., in spher- 
ical coordinates, V(r,9,Lp) V{r). The Hamiltonian 
is then S0(3) invariant, and its eigenfunctions, with 
good angular momentum quantum numbers, factorize 
into products fni{r)Yim{0,(fi) of radial wave functions 
and S0(3) spherical harmonics. Similarly, for the Bohr 
Hamiltonian, which is given in terms of quadrupole de- 
formation variables /3, 7, and (Sec. 2.2) by 
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(3) 



if the potential is a function V{(5) of the radial coordi- 
nate only, then the five-dimensional analog of a central 
force problem arises. The Hamiltonian is SO (5) invari- 
ant, and its eigenfunctions, with good seniority and angu- 
lar momentum quantum numbers, factorize into products 
fnv{l3)^vaLM{'y,^) of radial, i.e., /3, wave functions and 
S0(5) spherical harmonics, as in Ref. [18]. An example of 
application of the present methods to such "7-unstable" 
problems is found in Ref. [14]. 

A limited set of S0(5) D S0(3) spherical harmonics was 
computed many years ago by Bcs [19], for values of the 
angular momentum L < 6. These were obtained by se- 
ries solution of coupled differential equations in the coor- 
dinate 7, to find the eigenfunctions of the S0(5) Casimir 
invariant. However, this approach becomes prohibitively 
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complicated for L > 6. The S0(5) D S0(3) spher- 
ical harmonics arc the angular wave functions for the 
five-dimensional quadrupole harmonic oscillator [20, 21]. 
Thus, alternative approaches based on the oscillator ba- 
sis construction are possible, as summarized in Ref. [7]. 
A method based on the Cartan-Weyl reduction is given 
in Ref. [22]. 

Here we make use of the construction proposed in 
Ref. [10], in which the S0(5) D S0(3) spherical har- 
monics are developed as polynomials in a set of four 
basic generating functions defined on ^4. These func- 
tions were identified as generators of a complete linearly 
independent basis of SO(3)-coupled wave functions for 
L^{S4) [23], based on a knowledge of L'^{Si) as a di- 
rect sum of irreducible S0(3) subspaces [obtained from 
the S0(5) S0(3) branching rules]. The generating 
functions can be related to the "elementary permissi- 
ble diagrams" of the well-known algorithm of Chacon et 
al. [20, 21] for construction of a U(5) D S0(5) D S0(3) 
basis for the five-dimensional harmonic oscillator, imple- 
mented in the nuclear collective model code of Hess et 
al. [24, 25]. The salient property of the generating func- 
tions for L^(S'4) [10] is that they yield a direct route to 
the construction of S0(5) D S0(3) spherical harmonics 
without reference to the full harmonic oscillator problem. 
The method is set forth more concretely in Sec. 2.3. 



2.2. Representation of functions on S4, 

The quadrupole moments qm (m = 0, ±1, ±2) for the 
collective model transform as components of an L = 2 
spherical tensor under S0(3) rotations, i.e., 



(4) 



where ft represents the Euler angles for the S0(3) ro- 
tation, and the Wigner ^ function [26] is the rotation 
matrix element. [The quantities qm can alternatively be 
taken to represent the nuclear surface deformation pa- 
rameters ttm (rn = 0, ±1, ±2). The difference is only in 
physical interpretation and does not affect the following 
results for S0(5).] 

These quadrupole moments are conveniently ex- 
pressed in terms of Bohr's spherical polar coordinates 
(/?,7,0) [27], by the relation 



/3cos7^,i'i(l^) + ^/3sin7 



-2,m 



in) 



The squared length of a vector q e 



(5) 

is given by 



X^ml?™!^ = Thus, /? is the radial coordinate for 
and (7, fl) are angular coordinates. 

For consideration of the angular functions on the unit 
sphere S4, wc henceforth set /? = 1 and restrict considera- 
tion to the unit length quadrupole moments Qm, defined 



as 



Qr, 



cos7^£(f^) + -^sin7 



(6) 

These unit-length quadrupole moments are then propor- 
tional to the basic v = 1 spherical harmonics on ^4. 

Consider a function ^^'(7, Q) on S4, of good S0(3) 
angular momentum L and S0(2) quantum number M. 
Any such function may be expanded [27] in the form 



even 



(7) 



where 



{1 + 6kV/^ 



-KM\^'-) 



(8) 

Note that the 3l functions occur only in the symmetrized 



linear combinations 



'KM 



with even values 



{L) 
-KM 



KM 

{—)^s!'km, we can restrict 



of K . Thus, because ^ 

to > 0. Note also that the functions Ckm vanish 
identically for K = Q when L is odd. The normalization 
factor (1 + (5k)~^^^ (where 5k = 5k,o) is included for later 
convenience. 

The functions ^^2/ (^)> with K >Q, provide an orthog- 
onal basis for those fimctions of the S0(3) angles which 
respect the symmetry properties of collective model wave 
functions. Prom the inner product for the ^ func- 
tions [26], 



8^2 



2L-h 1 



-5l'l5k'k5m'm 



, (9) 



we obtain the inner product for the nonzero ^^lii^) 



£.^^?M'm'KU^)d^ = 



(L) 



2L + 1 



Sl'lSk'kSm'm- (10) 



The volume element on ^4 is given by dv = sin 37 dj dfl, 
where 7 is integrated over the range [0,7r/3]. Thus, the 
inner product (^'2 1 * 1 ) = / \E'2 * 1 sin 37 ^7 c/O of two func- 
tions on S4, when expanded according to (7), is given 
simply by 



16^2 



J2F^Kil)F2Kh) 



.k=o 

even 



sin 37 dj. 

(11) 



2L-I- 1. 

Functions of distinct L or M are orthogonal. 

2.3. SO(3)-coupled basis functions 



The generating function method [10] for construction 
of the S0(5) D S0(3) spherical harmonics rests upon the 



observation that a complete set of normalizable S0(3) 
highest- weight {i.e., those for which M = L) functions on 
Si is provided by the products 



*[ni,n2,n3,n4](7, 

= [$1(7, [$2(7, [*3(7, 



(12) 



where the exponents ni, 712, and 713 take on the vahies 0, 
1, . . ., and is restricted to or 1. These ^ [n i. ,12, na.rn] 
are monomials in four generating functions $1, $2; ^3, 
and $4, which are simply the M = L components 

$1 oc Q2 

*2 « (Q X Q)(') 



$3 oc (Q X Q X Q)^' 



(0) 



(13) 



$4 oc (Q X Q X Q)^' 



(3) 



of SO(3)-coupled products involving Q. The coupled 
product of two tensors is defined by 

yj XJ Jm — 2^ \ Ml M2\ M M2 -'■Ml ■ \^*) 

M1M2 

Note the right-to-left coupling order in this definition 
of the S0(3) tensor product, used for consistency with 
Refs. [8, 10] and to simplify the phases [28] arising in 
the Wigncr-Eckart theorem. The norms of these gener- 
ating functions are arbitrary and can be chosen for con- 
venience. In the standard form (7), we take 



$i(7,r!) = cos7 ^^^^(n)- 
$2(7, f^) = cos 27 ^^'^'(f^) - sm2j^g\n) 

$3(7,f^) = cos37C(f^) 

$4(7, ^i)= sin37d3^(^^). 



sin 7 ^g\n) 



(15) 



The multiplication of two highest-weight functions 
yields another highest-weight function (equivalent to the 
stretched coupling of two angular momentum functions). 
Hence, every $[„i,„2.n3,n4] ^ highest-weight function, 
with L = M = 2ni + 2n2 + 3^4 and "degree" N = 
n\ + 2n2 -|- 3n3 -|- 3n4 in the unit quadrupole moments 
Q [29]. The multiplication of the basis fimctions is car- 
ried out using the multiplication rule for the highest- 
weight functions S^xl- ^^^^ rivXe, which follows from the 
multiplication rule [26] for S> functions 



_\^( Li L2 
~ V Ki if 2 
L 



L W Li L2 I i V 
if A Ml M2\M }■ 



(16) 
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FIG. 1: The set of basis monomials $jvti of degree N <6, 
with labels A^, L, t, and K indicated. The numbered arrows 
indicate the raising actions of multiplication by the generat- 
ing functions $1, $2, $3, and $4. The angular momentum 
multiplicity at N = 6 is highlighted (dotted box) . This same 
diagram enumerates the "f^aL with v <6, that is, it gives the 
branching of S0(5) irreps {v, 0) into angular momenta L, if 
the axis label N is read as v. 



the equation 



AL2) 

'iK2L', 



SifiLi 



E 

if>0 



(l+<5if)l/2 



(1+(5kJ1/2(1+5^J1/2 



/ Li L2 \ Li+L2\j_(_\L2 ( Li L2 \Li+L2\ 
l>ifiif2| if ) l.ifi-if2l if ) 



fL\+L2 
?if,Li-|-L2- 



(17) 



where K = Ki-\- K2 and M = Mi -\- M2, is expressed by 



The more general tensor-coupled product of the spherical 
tensors which have components ^^fM) is considered 
in Sec. 2.4. 

An equivalent but more descriptive labeling for the 
highest- weight monomials $[„^,„2,„3,„^](7,0) is given by 
regarding each of them as the M = L component of the 
set of angular momentum L functions ^NtLMil,^) to 
which they naturally extend. Here monomials of the 
same degree TV and angular momentum L are distin- 
guished by the label t = n^, which counts the number 
of zero-coupled triplets of the quadrupole coordinates. 
In addition, the $ may be organized into quasibands of 
given bandhead angular momentum K ~ 2n2 + 2n4, as 
shown in Fig. 1. Together these functions ^NtLMil,^) 
constitute a spherical tensor $jvtL(7, ^)- Once the high- 
est weight function ^NtLLij,^) has been defined, the 
remaining ^NtLMil,^) follow immediately, since they 
share the same coefficients -Fif (7) in the expansion (7). 

An important characteristic of the generating function 
construction is that the set of L values appearing at a 
given degree N is identical to the set of L values arising 
for S0(3) irreps in an S0(5) irrep of seniority v = N. 
These are given by the known S0(5) D S0(3) branching 
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TABLE I: The SO(5) D SO(3) spherical harmonics '^vaL for v <3. For these lowest-seniority spherical harmonics, orthonor- 

malization is trivial, and the spherical harmonics are simply proportional to the corresponding monomial basis members 
{^vaL = ^ya'i^^'^^NtL)- The relations between the ^vaL (with L = and 2) and commonly-encountered scalar and quadrupole 
operators are given in brackets. 



^^r [=/§] 

41 cos 7 ^^^' + 4^ sin 7 [=^/fQ] 
^ (7 + 5 cos 27) 4 + 2^ sin 27 f ^ ^ + ^ ( 1 - cos 27) 
4Icos274f - 4Isin27ef [= -^(S x Q)*^)] 
A y/^(7 cos 7 + cos 37) ^^'^ + ^ ^(15 sin 7 + 11 sin 87) ^i^^ 

+ ^\/^(^°^ 1 - cos 37) if^ + ^ ^/f (sin 7 - sin 37) ^ 
I y^(5 cos 7 + 7 cos 37) C^'' - f sin 7 g*^' + | y^(- cos 7 + cos 37) 
I sin 37 

I cos 37 cos 37] 



V L ^vcl/{8it^) (87r2)i/2*„„i 

2/3 

1 2 4/15 

2 4 16/105 

2 4/15 

3 6 32/315 

4 88/945 

3 8/63 
4/9 



rule [16, 17]. Thus, Fig. 1 also enumerates the labels of 

all highest-weight spherical harmonics "^voLL np to the 
maximum seniority (v = 6) shown. The multiplicity cLnl 
of angular momentum L at degree N (also the multiplic- 
ity dyL of angular momentum L at seniority v) is given 
by (Al). A multiple occurrence of the same i at a given 
N first arises at A'' = 6, for L = 6 (Fig. 1). 

Although the ^NtLL form a complete set of highest- 
weight functions and are angular momentum eigenfunc- 
tions, they are, in general, non-orthogonal (for a given L) 
and also not eigenfunctions of A^. Thus, they are not the 
desired S0(5) spherical harmonics. Construction of the 
spherical harmonics 'i'vaLAiij,^) relies on the observa- 
tion that the highest- weight spherical harmonics '^voLL, 
like the ^NtLL, form a complete set of highest-weight 
functions and that they are polynomials of degree v in the 
unit quadrupole moments Q. For any given A^maxj the 

sets {"bvaLLlv < A^max} and {^NtLLW < A^max} both 

span the space of highest-weight polynomials of degree 
< Amax in Q. Construction therefore proceeds induc- 
tively, by orthonormalization of the bases of successively 
higher A^max- The 'i'vaLL, as polynomials of degree v in 
Q, are linear combinations of the ^NtLL with N < v. 
Because they must be orthogonal to all '^v'a'LL of lower 
seniority v' <v, they can therefore be obtained by Gram- 
Schmidt orthogonalizing the monomials of degree N = v 
with respect to the space of lower degree. 

The orthogonality of the spaces with differing L implies 
that orthogonalization can be performed separately in 
each space of given L. Within an i-space, the spherical 
harmonics up to seniority Wmax are obtained as follows: 

(1) Order the basis monomials ^NtLL, for < A^ < 



^^max, by increasing N, and then by increasing t when 
multiple basis monomials of the same L occur for a given 
A^ [30] . These functions may then be labeled with a single 
counting index, as ^u, with i = 1, . . ., D^^^^l- The 
dimension D^^^^l of the seniority-truncated L space is 
given by (A2) . 

(2) Calculate the overlaps for 1 < i,j < 
Dv^^^L, using (11). 

(3) Determine the linear transformation necessary to 
bring the into an orthonormal set by the Gram- 
Schmidt procedure. The result is a matrix of orthogo- 
nalization (and normalization) coefficients T^jj for the 
i-space, in terms of which the 

^Li= £ TLij^Lj (18) 

i=i 

are the desired spherical harmonics, in order of increasing 
seniority. 

In step (1), if angular momentum multiplicity occurs at 
a given A^, the ordering of the ^NtL sharing the same N 
is, in principle, arbitrary. However, choosing a different 
ordering for the $jvtL gives rise, after Gram-Schmidt or- 
thogonalization, to a different, equally valid, set of spher- 
ical harmonics at the corresponding seniority {v = N). 
That is, the spherical harmonics "^vaL (a = 1, . . ., dyL) 
span a seniority-degenerate subspace, and the ordering of 
the $7vtL in the orthonormalization process determines 
which of the possible unitarily-equivalent bases for the 
subspace is selected as the "spherical harmonics" . 

Note that the overlaps of the ^nill obey a parity se- 
lection rule. Since ^NtLL is a product of A^ factors of 
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Q, it has parity {—)^ under the parity operation, 
which takes — > —Qm- The overlap of two functions 
of opposite M'''-parity [31] vanislies, so {^N'L'LL\^NtLL) 
is nonzero only ii N + N' is even. The '^vaLL resulting 
from the orthonormalization process therefore retain def- 
inite R^-parity (— )" and are constructed only from the 
^NtLL of this same parity. 

Hence, by the W parity selection rule, all the ^NtLM 
with TV < 3 arc already orthogonal and therefore are 
the spherical harmonics with v < 3, to within normal- 
ization. The necessary normalization factors ^aL, such 

— 1/2 

that 'i'vaL = --^vaL ^NLL, and resulting spherical har- 
monics, for V <3, are summarized in Table I. The table 
also provides a glossary relating the '^vaL to Hamilto- 
nian (L = 0) and electric quadrupole (L = 2) operators 
commonly referenced in physical applications. More gen- 
erally, even for i; > 3, the two lowest-seniority ^'^ql for 
a given value of L are proportional to the ^NtL with 
N = v. Other higher-seniority "ifyaL are nontrivial linear 
combinations (18) of the $jvtL with N <v. 



2.4. Triple overlap integrals 

For the determination of S0(5) Clebsch-Gordan co- 
efficients (Sec. 2.5), we need to calculate triple overlap 
integrals 

(*3|*2|*i) = j *;(7,0)*2(7,fi)*i(7,0) sin37rf7rfO. 

(19) 

In this expression, ^'2 is interpreted as an operator which 
acts multiplicatively, i.e., *2|^i) has as its wave func- 
tion ^'2(7, fi)^'i(7, f2). For angular momentum coupled 



functions, we need only the reduced matrix elements 
(*^^=''||\E'^^'^||\E'^^'') defined by the Wigner-Eckart theo- 
rem [26] 



\*3M3l*2M2l*lMi/ 



(2L3 -I- 1)1/2 V Ml M2 I M3 ; 



X (* 



(i3)n^(i2)|i,T,(ii) 



(20) 



where the quantity in parentheses is an S0(3) Clebsch- 
Gordan coefficient. 

The Wigner-Eckart theorem is easily inverted {e.g., 
Ref. [8]) by application of the Clebsch-Gordan unitar- 
ity condition, to give the reduced matrix element in a 
computationally convenient form in terms of the coupled 

action of on Namely, 



{L3 



(i2) 



(i2) 



X I* 



{L3) 



M3 



(21) 



where [*^^'^ x has wave function 



,(Ll)vl(i3) 



*(^^)(7,l^)xvl/(^^)(7,l]) 



(is 



Ms 



E (M;M2|M3)*2Ml(7,0)M^^)(7,fi), (22) 

M1M2 

following the convention (14) for the coupled product of 
tensors. The functions $^m (-^ = 0, ±1 



±L), re- 



garded as components of a spherical tensor , obey the 
coupling rule 



,1/2 ( 



Li L2 



(1 + 5kJ V2(l + <5^Jl/2 \K,K2\ KM I^K,+K2 



+ 



{i+SK.-K.y/' j(-)^-(^^ ^ 

-5xJl/2(l+5^Jl/2 

obtained by direct application of (16). It follows that [33] 



Ki>K2 



K2 I K1-K2 )^Ki-K2(^) 
(l+5^jV2(l + 5^Jl/2 \i-)^^{Xk\-K'+K2)^%+K2m K,<K,^ 



(23) 



(i3)||^(i2)||jiy(il)\ _ 



') = 



167r2 f 

^3 + 1)1/2 J 



(2L3 + 1)1/2 



E 



even 

Li I Li L2 \L3\ Ki< K2 



(_\Li ( Li L2 
^\ J \ -Ki K2 



K3, 



-Flifi (7)^2 1^2(7)^3X3 (7) 



sin 37^7. (24) 



2.5. SO(5) D SO(3) Clebsch-Gordan coefficients If the S0(5) representations are labeled according to the 

S0(5) D S0(3) D S0(2) subalgebra chain (1), each cou- 
Representations of S0(5) couple to form new represen- pling coefficient may be written as the product of an 
tations according to S0(5) Clebsch-Gordan coeflacients. 
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S0(3)-rcduccd Clebsch-Gordan coefficient and an ordi- 
nary S0(3) Clebsch-Gordan coefficient, according to the 
Racah factorization lemma [34] (see Ref. [35] for a gen- 
eral discussion). We term the S0(3)-reduced factor an 
S0(5) D S0(3) Clebsch-Gordan coefficient. Couphng co- 
efficients, such as these, which are reduced with respect 
to a subalgebra are also known as "isoscalar factors" [36] . 

In the context of the S0(5) spherical harmonics, 
only the symmetric representations {v, 0) of S0(5) arise. 



These are related by the S0(5) D S0(3) Clebsch-Gordan 
coefficients ( ^"'f^ ^"^f^ I ) , where the ai, aa, and as 
are multiplicity indices, corresponding to the multiplic- 
ity in the Let l°mJ and {Xa^L%J denote 
orthonormal bases for S0(5) representations of seniority 
Vi and V2, respectively. Then an orthonormal basis for 
an S0(5) tensor-coupled product irrep of seniority v is 
defined by 



ai ,Li,a2^L2 
Mi,M2 



(^^1,0) (^^2,0) 

aiLi aa-Zva 



iv,0)\ 




i2 




aL J 




Ms 


M J 



(''2,0) 
Aa2i2M2 



(-1,0) 
AaiLiMi 



(25) 



The S0(5) D S0(3) Clebsch-Gordan coefficients en- 
ter into the S0(5) Wigncr-Eckart theorem, which gov- 
erns matrix elements of S0(5) D S0(3) tensor operators. 
The Racah factorization lemma may be used to write the 
S0(5) Wigner-Eckart theorem in terms of S0(3)-reduced 
quantities as 



v^asL^ ||^f2a2l'2 ll^«ic«i-£/i) 

(^^1,0) {V2,Q) 



(^^3,0) 

asLs 



(26) 



where (*„3 [j^^^ is an S0(5)-reduced_ (doubly- 
reduced) matrix element. The factor \/2Li + 1 compen- 
sates for the corresponding factor absorbed into the defi- 
nition of the S0(3)-reduced matrix element in (20). The 
S0(5) D S0(3) Clebsch-Gordan coefficients satisfy the 
normalization condition 



E 

aiLi 



(1-1,0) (^;2,0) 
aiL-i aa^a 



(1^3, 0) 

as 1/3 



= 1 



(27) 



from unitarity, where the phases of the Clebsch-Gordan 
coefficients are real by convention. The values for L 
occurring within a given S0(5) representation (fjO), 
and their multiplicities, are governed by the multiplic- 
ity formula (Al). The Clebsch-Gordan coefficients van- 



ish unless Li, L2, and Z/3 satisfy the triangle inequal- 
ity {\Li — < Ls < Li + L2). The values of vi, 
V2, and V3 likewise must satisfy the triangle inequality 
— ^2 1 < ws < ui -|- W2), as well as the parity constraint 
that vi +V2 + V3 be even. 

The S0(5) D S0(3) Clebsch-Gordan coeffi- 
a sis' ) follows immediately from 
the corresponding S0(3)-rediiced matrix element 
{'^v3a3L3\\'^v2a2L2\\'^viaiLi) , by (26), once the S0(5)- 
reduced matrix element (^'„3||4'„2||^'„i) is known. In 
fact, from (26) and the normalization condition (27), we 
obtain 



cient ((--o)("-o) 

V aiLi OL2L2 



(^--3 11*^2 ll*-l)^= E 



I* 



aiLi 



2Ls + l 



(28) 

for each L3 and a^. This yields the S0(5)-reduced matrix 
element in terms of the summed squares of the computed 
S0(3)-reduced matrix elements, to within an arbitrary 
phase, which is chosen to be unity. 

Explicit closed-form expressions for the S0(5)- 
reduced matrix element have been given for the case 
(\I'„_l_i||\l/i||\I'„), i.e., for the quadrupole tensor [3]. More- 
over, by numerical inspection of the normalization sums 
for an extensive set of computed S0(3)-reduced matrix 
elements, we find [37] that the S0(5)-reduced matrix el- 
ements are given by 



II* II* )- 1 ^/ (2^i+3)(2.;T^ + 

\ .3 II V2\\ V,) (^3 + 2)(t;3 + l) {la - v,)l{^a - V2)Kla - vs)\ 



L^^,^i^-2v. + m^-2v2 + m^-2vs + iy.^ (29) 



(a + 3)! 
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with a= (v\ + U2 + vj^. This conjectured expression has 
been verified exhaustively for all combinations of values 
for wi, t;2, and U3 likely to be considered in nuclear col- 
lective model calculations. However, we stress that is 
has not been proved in full generality. Use of the ex- 
pression (29) for (*i,3 11*1,2 11*^^1) in tlic S0(5) Wigner- 
Eckart theorem (26) allows each S0(5) D S0(3) Clebsch- 
Gordan coefficients to be extracted directly from the cor- 
responding S0(3)-rcducc(l matrix elements (24), without 
the need to calculate matrix elements for all L\ and L2 
values involved in the normalization condition (27). [The 
use of (29) in the code is described in Sec. 4.2.] The va- 
lidity of the conjecture, for given V2, and v^, may 
be established by calculating the normalization sum for 
the imdcrlying computed S0(3)-reduced matrix elements 
or, equivalently, by explicitly verifying that the extracted 
coefficients satisfy (27). 

The S0(5) D S0(3) Clebsch-Gordan coefficients obey 
symmetry relations [10] 



(«2,0) («i,0) 



(«3,0) 



(-) 



L1+L2— i^a 



(^^1,0) (^^2,0) 

aiLi a2L2 



(30) 



and 



iv3,0) {V2,0) 

asL^, a2L2 



(«i,0) 
aiLi 



L1+L2 



d^,{2L3 + l) f {v,,0) {v2,0) 



dvs{2Li + 1) V aiLi a2-t'2 



(^^3,0) 

0:3 L3 



(31) 



where dv = ^{v + l){v + 2){2v + 3) is the dimension of 
the S0(5) representation {v,0), obtained from the Weyl 
formula {e.g., Ref. [35]). Hence, it is only necessary to 
calculate the Clebsch-Gordan coefficient for one permu- 
tation of {Vi ,V2,Vs). 



3. IMPLEMENTATION 
3.1. Overview 

The computer code for construction of S0(5) D S0(3) 
spherical harmonics and calculation of S0(5) D S0(3) 
Clebsch-Gordan coefficients is implemented set of 
packages for Mathematica 6 [15]. Mathematica provides 
native support for symbolic arithmetic. In the present 
context, this allows all calculations to be carried out ex- 
actly, in terms of expressions involving square roots of 
rational numbers. 

The basic algorithm for the calculation of S0(5) D 
S0(3) Clebsch-Gordan coefficients, as described in Sec. 2, 
involves three main computational tasks: 

(1) construction of the monomials (12) comprising the 
basis of highest- weight functions, 



(2) calculation of the overlaps (11) of these monomials, 
as needed for the orthonormalization process, and 

(3) calculation of triple overlaps (24), as needed for the 
S0(5) D S0(3) Clebsch-Gordan coefficients. 

For practical implementation, two algorithmic refine- 
ments are made to this scheme. These relate to the in- 
ternal representation of the fimctions Fk{^i) (Sec. 3.2) 
and to the recognition (and full utilization) of extensive 
redundancies among the calculations involved in evaluat- 
ing different triple overlap integrals (Sec. 3.3). Together, 
an efficient treatment of these aspects of the implementa- 
tion extends the range of applicability of the method (for 
calculation in exact arithmetic) from a maximal seniority 
of approximately 10 to seniorities of > 100, thereby eas- 
ily yielding more than sufficient S0(5) D S0(3) Clebsch- 
Gordan coefficients for nuclear structure calculations. 

To understand the considerations underlying the effi- 
cient implementation of the algorithm, let us more closely 
consider the structure of the calculation. The computa- 
tionally most demanding task is evaluation of the over- 
lap or triple overlap integrals in the '^nil basis. Re- 
call that this involves first constructing the function of 
7 appearing in the integrand of (11) or (24) and then 
evaluating the integral with respect to 7. The inte- 
grand is built from the generating functions (15), us- 
ing (12) and (17). The integrand is thus a polynomial 
in the trigonometric functions cos 7, sin 7, cos 27, sin 27, 
cos 37, and sin 37 or, therefore, by multiple-angle iden- 
tities, a polynomial in cos 7 and sin 7. For the overlap 
{^N2t2LM\^NitiLM), the resulting polynomial requires 
powers of cos 7 and sin 7 as high as N\ + N2, or, for 
the triple overlap {'^N:iL?.L?.\\^N2t2L2\\'^NitiLi) , as high 
as Ni + N2 + N3. The integral may be evaluated ex- 
actly, but simplification of the polynomial and integra- 
tion with standard symbolic algebra software becomes 
prohibitively inefficient for the construction of spherical 
harmonics with v > 10. Alternatively, attempts at brute- 
force numerical integration are hampered by the highly- 
oscillatory nature of such polynomials in trigonometric 
functions. Thus, it is seen that evaluation of the overlap 
integrals is the defining computational challenge. 

An effective solution arises from the realization that 
any polynomial of degree n in cos 7 and sin 7 can be de- 
composed as a finite sum of exponentials 



/(7) = E 



Ik-y 



(32) 



k=—n 



i.e., as a finite Fourier series, of degree n. Integration of 
exponentials is trivial. Such a Fourier expansion method 
was, in fact, used in the final version of the code used 
to calculate the Clebsch-Gordan coefficients tabulated 
in Ref. [10]. In the present implementation, from the 
very beginning of the calculation, we simply represent 
all ^V(7) coefficients, starting with those of the integrity 
basis functions (15), as finite Fourier series (32). That 
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is, every function is replaced by a list of Fourier coeffi- 
cients {a_„, . . . , flg, . . . , o„}. Then, when the integrand 
in (11) or (24) is constructed as a function of 7, it is al- 
ready manifestly Fourier expanded, and the 7 integration 
is straightforward. 



3.2. Fourier series representation of functions 

When the functions of 7 are represented as Fourier 
sums (32), there are three basic arithmetic operations 
which must be carried out on the corresponding sets 
of Fourier coefficients. Addition of two functions f{j) 
and 5(7) is accomplished by addition of their Fourier 
coefficients. Multiplication of a function by a constant 
is simply accomplished by multiplication of the coeffi- 
cients by that constant. Multiplication of two functions 
7(7) and 5(7) by each other gives rise to a convolu- 
tion of the coefficients. Specifically, let 7(7) and (7(7) 
be given by Fourier sums 7(7) = J2T=-m^r{e^^Y and 
5(7) = J2^=-n ^s{e-^'^Y ■ Expanding the product and col- 
lecting like powers of e^'^ yields the product rule 



/(7)5(7) = E '^'^(e^')' 



(33) 



fe= — (m+n) 

with the new Fourier coefficients given by 

min(2m— /c,2n+fc) 
t=max(— 2m— fe,— 2n+/c) 



(34) 



This is essentially the Fourier convolution theorem, in 
discrete form. 

For real-valued functions which are even in 7, the co- 
efficients Ofe are pure real and obey the symmetry con- 
dition a_fe = flfe. Similarly, for odd real- valued func- 
tions, the coefficients are pure imaginary and obey the 
symmetry a_fc = — a^. For instance, the functions aris- 
ing in the definition of the generating functions (15) are 
C0STO7 = [(e+'T)" + (e-'T)'"]/2 and sinTO7 = [(e+'T')" - 
(e~*''')™]/(2i), which arc even and odd, respectively. In 
either case, only coefficients with fc > need be stored, 
and only real-valued coefficients are required, provided a 
factor of i is absorbed into the definition of the series for 
odd functions. 

Thus, we use the representation 



(1 \ ff " 



(35) 



with 3 = for even functions and g=l for odd functions. 
This is effectively a Fourier cosine series or a Fourier sine 
series, for the even and odd cases, respectively. 

The definite integrals needed for the calculation are of 
the form j^^^ f(pi) dq. Moreover, only odd integrands 



[g = 1) arise in the problem. The definite integral of an 
odd series (35) on the "sector" 0<7<7r/3 is 



/■7r/3 

io • 



/(7)rf7 = E 



fe=i 



2ak 



1 — cos 



{k mod 6)77 



. (36) 



Evaluation of this sum requires only a limited set of 
trigonometric values, namely, cos kn/S {k = 0, 1, . . ., 5). 

The necessary definitions for working with Fourier rep- 
resentations of functions are contained in the subpackage 
FourierSum. The Fourier sum (35) is represented sym- 
bolically by the expression 

FourierSum [BohrGainma,5f,{ao,Oi , . . .,o„}] , 

where g indicates the symmetry {g = or 1), and the Uk 
are the real- valued Fourier coefficients as defined in (35). 
(The tag BohrGamma simply serves to indicate that the 
argument is 7, since the package allows for more general 
possibilities.) Thus, for example, cos 27 is represented as 

FourierSum [BohrGamma, 0, {0,0, 1/2}] . 

The coefficients which arise in the present calculations 
are rational numbers or square roots of rational numbers, 
and they are maintained as exact symbolic expressions 
throughout the calculation. 

The package defines addition of two FourierSum ex- 
pressions, multiplication by a constant, and multipli- 
cation of two FourierSum expressions as extensions to 
the usual Mathematica + and * operations. The func- 
tion IntegrateSector [r ,/] returns the integral given 
in (36). The function FourierSumToTrig [/] is also pro- 
vided, to convert FourierSum expressions back into ordi- 
nary symbolic expressions in terms of trigonometric func- 
tions. 



3.3. Evaluation of matrix elements 

The task of computing matrix elements (triple over- 
lap integrals) of the S0(5) D S0(3) spherical harmonics 
is most conveniently carried out by first calculating the 
matrix elements of the original monomial basis functions 
^NtL, that is, the ($^3*3^3 11*^2^*2^2 II *Witiii)- The de- 
sired matrix elements {^vsasLa \\'^V2a2L2\\'^viaiLi) in the 
orthonormal spherical harmonic basis then follow imme- 
diately from the Gram-Schmidt transformation (18). It 
is simplest to write this using the counting index labeling 
within an L-space, as 

(4-^3,3 ||*L2,2ll*Lm) = 

E '^LatijiTL^i^h^Lthni^LajsW^LihW^Liji)- (37) 

3l32j3 

As described in the preceding section, we have trans- 
formed the relatively intractable integration problem into 
the more manageable task of constructing products of 
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Fk{i) functions in Fourier representation. Therefore, 
the computational burden now hes primarily in evaluat- 
ing the discrete Fourier convolutions (33 34). Although 
the operations involved are in principle simple arithmetic, 
the number of iterations required for a single convolu- 
tion is substantial, growing as the square of the degrees 
of the Fourier sums involved, and the arithmetic itself 
involves symbolic simplification of expressions involving 
square roots of rational numbers. 

The challenge, therefore, lies in minimizing the number 
and complexity of the Fourier convolutions which must 
be evaluated. The problem is best approached by noting 
that the matrix elements are not to be calculated singly, 
but rather in aggregate, that is, as the set of all matrix 
elements of an operator ^N2t2L2 between i-spaces Li and 
L^. The integrand in the expression (24) for the matrix 
element is a sum over products of ^^(7) functions, and 
products of the exact same pairs of Fk{j) functions arise, 
redundantly, in the evaluation of many different matrix 
elements. 

A few straightforward observations allow us to remove 
these calculational redundancies. Most obviously, ev- 
ery matrix element (^^WataLall^WstsLall'I^WitiLi) involv- 
ing the same ket |$jVitiLi) and operator ^i^^t^^^ will 
also involve the same products between i^x (7) functions 
from '^NitiLiil-,^) and $7V2t2L2 (7i the notation 

of (21), it is therefore advantageous to calculate the cou- 
pled action of ^N^t^L^ on any given |$jvitii,i), namely, 
[$^2(2^2 X I'&AriiiLi)]^"^-''-'; only once, and to reuse this 
intermediate result for the matrix element with each bra 

More significant, though, is the observation that all 
powers of $3 cx cos 37 factor out of the summations 
in (21). (This will be more clearly apparent below.) 
Thus, the integrands involved in the calculation of many 
different ($Ar3t3L3||^'Ar2t2L2ll*AfitiLi). sharing the same 
total label t = t\ + t2 + tz, are actually identical. Even 
greater reduction in the number of convolutions needed 
is obtained by first evaluating the integrand only for 
those monomials involving no powers of $3 (ti = t2 = 
tz — 0) and only then multiplying by the relevant power 
(cos 37)*. 

In practice, the entire process is based on operations 
involving the coefficient functions FKiri) which describe 
functions on ^4, as in (7). The relevant definitions are 
provided by the subpackage WignerDSum. A function on 
^4 is represented as the expression 

WignerDSum [i,{Fo,i^2,...,i^[Lj 2}] > 

where L is the angular momentum and \ L\ 2 denotes the 
greatest even integer less than or equal to L. The in 
turn, are FourierSum expressions (Sec. 3.2). Thus, for 
example, the generating function $1 from (15) is repre- 
sented by 

WignerDSum [2 , { 

FourierSum [BohrGamma.O, {0,1/2}] , 
FourierSum [BohrGamma, 1 , {O , 1/2}] 

}] • 



Operations of addition of two functions and multiplica- 
tion by a constant are readily defined, by entry-wise op- 
erations on the lists of Fxici) coefficients. These oper- 
ations are defined by the package as extensions to the 
usual Mathematica + and * operations. 



To obtain the coupled action ^t^^t^h^ ^ 1^ 



l(i3) 



we evaluate the coupled product [$jv2t2i-2(7) ^) ^ 
^JViiiLi(7, ^^)]^^^^- The calculation follows in a 
straightforward fashion from the coupling rule (23) 
for the functions Cjc ^- The coupling operation is 
provided by the WignerDSum package as the function 
TensorCouple , ^2 , • 

However, the monomials $jv2t2/^2(7) ^) and 
^JVitiLi(7)f^) niust first themselves be constructed, 
from the definition (12). As noted in Sec. 2.3, mul- 
tiplication of highest-weight functions is identical to 
stretched {L = Li + L2) tensor coupling. The stretched 
product rule (17) is simply the L = Li + L2 special 
case of the more general coupling rule (23) for the 
Therefore, the $jvtL are also constructed in the code by 
use of this same function TensorCouple, as 



NtL 



(2)1 



(3)] "4 
4 J 



(2) 

where is the tensor with highest-weight component 
$1(7, n), etc., and L = 2ni-|-2n2+3n4 (Sec. 2.3). To com- 
plete the calculation of the matrix element, the overlap 
with {^NstaLsl must be evaluated. In fact, the function 
TensorCouple also suffices for this purpose. The sum 
of products of Fk{'^) arising in the overlap integral (11) 
is readily expressed in coupled form. By the coupling 
rule (23), 

even 

(39) 

This follows from the basic expression for L = 

Clcbsch-Gordan coefficients, much like the usual 

scalar product result [t!^^^ x t[^^\^q^ = (-)^(2L + 

1\-1/2Y^ ( \M rp(L) rp(L) 

It remains then to carry out the integration over 7, by 
use of IntegrateSector (Sec. 3.2). If we define {f{j)) = 

^''^'^ lo^^ fil) sinS-yd-y, where the factor of Stt^ results 
from the integration over Euler angles, then the reduced 
matrix element deduced from the overlap (21) is simply 

{^NatsLs II *JV2t2i2 II ^JVitiLi ) 

= ([^iVstsLs X {^N2t2L2 X ^iV.t.Lj'^^^^lD- (40) 

Expressing the overlap in terms of a coupled prod- 
uct docs not change the underlying calculation (and at 
most affords a convenient economy of coding). However, 
this formulation does permit the factorization of (cos 37)* 
from the integrand to be expressed in an especially sym- 
metric (and explicit) form. Observe that $3°^ oc cos 37 is 



(38) 
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a scalar, and, therefore, unlike the other generating func- 
tions, it factors out of the coupled product defining the 
monomial in (38), and subsequently out of the couplings 
in the zero-coupled product in (40). Hence, if we label 
each of the monomials by its exponents with respect to 
each of the generating functions, the integral factorizes 



as 



[i/i, 1^2, 1^3, 1^41 [ni,n2,n3,ni 



^(i3)j(0) 



X ($ 



(i2) 

[1^1,1/2, 0,1/4] 



(ii) 

[ni, 712,0,714]^ 



AL3 



'](0) 
J 



X $ 



[0,0,t,0]) 



(41) 



where t 



n's + t^3 + ns is the total exponent of $3°' 



occurring in the product. Note that only monomials of 



the form $ 



ni ,n2 ,0,n4 



, with Us = 0, are actually needed in 



the calculation, along with $[o,o,t,o] = (v^ cos 87)*. 

The mechanics of the optimized calculation of matrix 
elements therefore proceed as follows. The monomials 
^[ni,n2,o,n4] ^re Constructed by (38). To eliminate re- 
dundancy, this is done recursively, that is, by multipli- 
cation of the monomial of one lower degree in rii, n2, 
or 714 by the appropriate generating function, and all 
the intermediate results are cached. Then the succes- 
sive coupled products in (41) are constructed and cached: 

($(^^) X , O^"^'^ (stage I) 

(*£l2,o,.4] X ^ln:Lo,n/"'V (Stage II), and finally 
the full integrand with $[0,0, t,o] included (stage III). 

For matrix elements within an L-space, Hermiticity re- 
duces the number of independent matrix elements nearly 
by half. Since, under complex conjugation, ^NtL satisfies 



[#'^ n '1 X 



NtLM 



{-)^-^^NtL-M, it follows [26] that 



{^NatsLs \\^N2t2L2 \\^NitT_Li ) 

= (-)^^+''^-'^M*^.t.Ljl^iV2t2L2ll<i'iV3t3L3)*- (42) 

As a specific example of the computational savings pro- 
vided by the caching of intermediate products, consider 

the calculation of matrix elements of 3^ii2 Q within 
the space Li = L3 = 40, for all monomials of degree 
< 50. [This calculation yields the S0(5) D S0(3) Clebsch- 
Gordan coefficients with V2 = 1 and £2 = 2, for all vi and 
< 50.] There are i?5o,4o = 154 basis functions and 
thus 23 716 matrix elements to be calculated, or 11935 
after the relation (42) is taken into account. Evahiation 
of the coupled action of <I>ii2 (stage I) involves only 21 
couplings, since the coupling is only carried out for basis 
functions with ris = 0. Evaluation of the scalar cou- 
plings (stage II) constitutes the bulk of the calculation. 
Whereas each of the couplings evaluated in stage I in- 
volves a low-L, low-degree factor ($112), in stage II both 



factors, $ 



(-^^3) 



and ($ 



(i2) 



X $ 



AL3 



[ni,n'2,0,n'^\ V"^ [1^1,1^2, 0,1/4] ^ "^[ni,n2,0,n4]/ ' 

are typically of high L (hence each contains many Fk 

terms) and high degree with respect to 7 (hence many 
terms are involved in the Fourier sums). A total of 231 



200 



eioo 



H 2 




FIG. 2: Basis sizes and corresponding execution times, for cal- 
culating the matrix elements of ^'112 oc Q between the spheri- 
cal harmonics with v < 50. (a) The number of basis functions 
^fmaxi (i^max ~ 50), foi L cvon. (b) Time for evaluating 
the matrix elements within a single L-space (Li — L3 = L), 
again only shown for L even. (The basis dimensions and, 
consequently, execution times are significantly smaller for L 
odd.) All times are for calculations carried out under Math- 
ematica 6 for Linux, running on a 2.2 GHz Advanced IVIicro 
Devices Opteron processor. 



such couplings are necessary. Many more distinct expres- 
sions, namely, 3001, must be evaluated at stage III, but 
each entails only a single multiplication of scalar func- 
tions (Fourier convolution) and is therefore not compu- 
tationally intc^iisive. 

Total execution times for evaluation of these matrix 
elements within each i-space, as well as the dimensions 
of these spaces, are shown in Fig. 2. It is seen that calcu- 
lations for higher L are essentially more demanding than 
those for lower L, even if the bases for the L-spaces are of 
the same size. For higher L, a larger number of Fk terms 
are involved in each calculation. Furthermore, the basis 
monomials tend to involve higher powers of those gen- 
crating functions which carry angular momentum {^\ , 



$2^'' , and ) , in preference to , lessening the extent 
to which the optimizations discussed above can simplify 
the calculation. 

The full, optimized process of evaluating the ma- 
trix elements of a monomial ^iVjtjLa between two 
L-spaces Li and L3 is carried out by the function 
CalculateMonomialMatrix. The considerations just de- 
scribed also apply to eliminating redundancies in the 
evaluation of the overlaps needed for the orthonormal- 
ization process. 



S(3) 



(0) 
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4. USE OF THE COMPUTER CODE 
4.1. Installation 

The Computer Physics Communications Program Li- 
brary deposit associated with this article contains Math- 
ematica package files GammaHarmonic.m, FourierSum.m, 
WignerDSum . m, GramSchmidt .m, and Cache. m, in ASCII 
text format. Annotated code is included for the pack- 
ages in both Mathematica notebook format and Portable 
Document Format. The deposit also includes a Math- 
ematica notebook demonstrating the use of the code 
(Excunple .nb) and example output files containing tabu- 
lations of Clebsch-Gordan coefficients (see Sec. 4.2). 

All package files must be placed in a directory in 
the Mathematica file search path, as explained in the 
Mathematica documentation [15], or else in the cur- 
rent working directory for Mathematica. The pack- 
age GammaHarmonic must then be loaded, by eval- 
uating Get ["GammaHarmonic"'] , as demonstrated in 
Example . nb. 



4.2. Instructions 

The GammaHarmonic package provides calculation con- 
trol functions to carry out the several steps involved 
in constructing and tabulating a set of S0(5) D S0(3) 
Clebsch-Gordan coefficients. The syntax for each of the 
control functions is given in Table 11. Here, we provide 
basic instructions for carrying out these calculational 
tasks. A worked example, with sample input and out- 
put, is given in the Mathematica notebook Example. nb. 

(1) The labeling scheme for the basis monomials must 
first be set up, using ConstructCjNSet (see Table II for 
the appropriate syntax). This function assembles and 
stores a list of (iV, L,t,K) labels for the monomials {N < 
Vma.^) spanning each L-space {L< Z/max) [38]. 

(2) The necessary monomials ^NtL for the calculation 
must be constructed, with ConstructMonomials. Only 
monomials with na = are evaluated, for the reasons 
discussed in Sec. 3.3. 

(3) The function ConstructMonomialMatrices is then 

used to construct matrix elements (^LgiaH ^1,2*2 II ^iin) 
in the monomial basis, by the process discussed in 
Sec. 3.3. The function calculates all matrix elements for 
a given L2 and ^2, between all L-spaces Li and L3, up 
to I/max, subject to the triangle inequality. This calcu- 
lation should be carried out for any (^2,^2) which con- 
tribute, in (37), to the final spherical harmonic matrix 
elements of interest. For nuclear structure applications, 
these would typically include (£2,^2) = (2, 1), for the elec- 
tric quadrupole operator, and (-^2,22) = (0,2), needed 
for cos" 37 contributions to the potential in the Hamil- 
tonian. The calculation must also be carried out for 
{L2, ^2) = (0, 1), i.e., for the identity operator <I>oiOj since 
this provides the overlaps needed for the Gram-Schmidt 



orthonormalization. 

(4) The function ConstructBasisGST must be used to 
generate the Gram-Schmidt transformation coefficients 

(5) Finally, the function ConstructBasisMatrices 

is used to calculate the reduced matrix elements 
(^igjg 11*^282 II ^i-iii) of the spherical harmonics, by 
Gram-Schmidt transformation (37) of the previously- 
calculated monomial matrix elements. 

The S0(5) D S0(3) Clebsch-Gordan coefficients fol- 
low from these computed S0(3)-reduced matrix el- 
ements, by the S0(5) Wigncr-Eckart theorem (26) 
and the normalization condition (27), as described in 
Sec. 2.5. Individual Clebsch-Gordan coefficients may 
be accessed with the function SOSClebschGordanCumax, 
{vi,Li,ai} ,{v2,L2,a2} ,{v3,L3,a3}'\ . So that the co- 
efficients can be extracted independently, without requir- 
ing computation of matrix elements for all Li and L2 
values encountered in the normalization sum (27), the 
function SOSClebschGordan directly applies the expres- 
sion (29) for the S0(5)-reduced matrix element in the 
S0(5) Wigner-Eckart theorem. 

The function WriteS05CGTable outputs a tabulation 
of S0(5) D S0(3) Clebsch-Gordan coefficients sharing the 
same V2, 0.2, and L2. These correspond to spherical har- 
monic matrix elements sharing the same ^«2a2i2- Coef- 
ficients are listed in order of increasing L3, Li, ^3, a^, 
v\, and ai (that is, with the last of these indices vary- 
ing most rapidly). Only coefficients with Li < L3 are 
tabulated, in recognition of the symmetry relation (31), 
and only those coefficients for which nonzero values arc 
allowed under the S0(3) triangle inequality and S0(5) 
selection rules (Sec. 2.5) are included. 

The Clebsch-Gordan coefficients arc written both ex- 
actly and as floating point numbers. The exact value 
of a Clebsch-Gordan coefficient can be expres sed as the 
signed square root of a rational number, ±^/a/b. The 
magnitude of each value is squared in the output, to elim- 
inate the need for radicals, so the value written is zLa/b. 
Because the numerators and denominators appearing in 
these rational numbers can be large (exceeding 400 deci- 
mal digits, for instance, for Umax = 50), the floating point 
value is also given. This is meant to facilitate input by 
external programs written in languages which do not pro- 
vide native support for arbitrary-length integers. Each 
row of the tabulation has the form 

vi Li ai V2 L2 Q!2 V3 L3 as X ± a/b, 

where x is the floating point representation of ±a/b. 

Tabulated coefficients may be read back with 
ReadS05CGTable. The corresponding S0(3)-reduced ma- 
trix elements are recovered and stored as matrices Basis- 
0peratorMatrix[i;max,{-^2.i2}.{i3.ii}], which may 
then be used, for instance, for ACM calculations in Math- 
ematica. When using the matrix elements calculated by 
the code, it should be noted that the factor of Stt^ which 
arises from the integration over Euler angles has been 
suppressed in the actual calculations. (This choice is 
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TABLE II: Control functions for S0(5) D S0(3) spherical harmonic and S0(5) D SO(3) Clebsch-Gordan coefficient generation. 



Function Description 



ConstructQNSet [Vmax .-C-max] 

ConstructMonomials [Vmax . imax] 
ConstructMonomialMatricesCvmax .imax.{i2 .12}] 

ConstructBasisGST [Umax . imax] 
ConstructBasisMatrices [Vmax .i/max .{i^2 .12}] 

Wr iteS05CGTable [Vmax , imax , {^2 , 42 } .filename] 
ReadS05CGTable [Wmax , imax , {^2 , 12 } .filename'] 



made in order to eliminate the symbol tt from input and 
output.) If the true normalization is desired, all com- 
puted matrix elements must be multiplied by (Stt^)"^/^. 

Two example tabulations of S0(5) D S0(3) Clebsch- 
Gordan coefficients are included with the Computer 
Physics Communications Program Library deposit. 
These contain coefficients with {v2, L2) = (1, 2) and (3, 0) 
[that is, (i2,*2) = (2,1) and (0,2)], for all seniorities vi 
and V3, up to Wmax = 50 and for all allowed values of Li 
and I/3 at these seniorities (Lmax = 100). The tabula- 
tions are given in the files basis-50-100-cg-2-l.dat 
and basis-50-100-cg-0-2.dat. (The numerical labels 
indicate fmax, -f'max, L2, and 12, respectively.) The tabu- 
lated coefficients are of sufficient extent to support basic 
nuclear structure calculations with the ACM. 

Intermediate results may be saved to files at various 
stages of the calculation, and later retrieved, allowing 
computations to be resumed or extended [39] at a later 
time. In particular, monomial matrix elements may be 
saved with WriteMonomialMatrices and subsequently 
read back with ReadMonomialMatrices (see Exsmiple .nb 
and the internal program documentation). Similar func- 
tions are defined to write and read back the Gram- 
Schmidt transformation coefficients and to write out the 
S0(3)-reduced matrix elements of the spherical harmon- 
ics in matrix format. 



4.3. Explicit expressions for the SO(5) D SO(3) 
spherical harmonics 

Explicit expressions for the S0(5) D S0(3) spherical 
harmonics, as functions of 7 and the Euler angles (Ta- 
ble 1), are not needed for diagonalization of the collective 
model Hamiltonian or calculation of transition matrix 
elements in the ACM framework. The S0(5) D S0(3) 
Clebsch-Gordan coefficients suffice. However, if ACM 
calculations are carried out in an SU(1, 1) x S0(5) basis, 
probability distributions P(/3, 7) for the wave functions 



Constructs a list of monomial {N, L, t, K) labels for each 
L-space. 

Precaches the $jvti, constituting the monomial basis. 

Calculates the matrix elements ('I'Lsia ||$L2i2 ll'S'^in ) ™ ^ho 
monomial basis, by evaluation of the triple overlap integrals. 
Constructs the Gram-Schmidt transformation coefficients. 

Calculates the matrix elements (^'£313 ||'I'i,2i2 ll^^in) the 
spherical harmonic basis, by Gram-Schmidt transformation of 
the monomial matrix elements. 

Writes a tabulation of S0(5) D SO(3) Clebsch-Gordan 
coefHcicnts, as defined in tlic text, including both the exact 
rational and floating-point values. 

Reads back a tabulation of S0(5) D SO(3) Clebsch-Gordan 
coefficients, as defined in the text. 



in coordinate space can then be obtained by combining 
the known radial wave functions -R^(/3) [3] and the ex- 
pressions for the ^'„ql(7, fi). The general procedure is 
given in the appendix of Ref . [9] , and example probabil- 
ity distributions may be found in Figs. 3 and 4 of that 
reference. 

The function ConstructBasisWaveFunctions explic- 
itly constructs the spherical harmonics '^vahil,^), by 
taking linear combinations of the monomials ^mL, ac- 
cording to the Gram-Schmidt transformation (18). The 
spherical harmonics may be converted from Fourier rep- 
resentation to symbolic expressions involving trigonomet- 
ric functions by use of FourierSumToTrig (Sec. 3.2), as 
demonstrated in Example . nb. A normalization factor of 
(Stt^)"^/^ must be restored in these results, as indicated 
in Table I, if the true normalization is desired. 



4.4. Algebraic collective model infrastructure 

If ACM calculations are to be carried out within 
Mathematica, the GciinmaHarmonic package provides sev- 
eral additional functions which can facilitate this pro- 
cess. As already noted (Sec. 4.2), the package provides 
the necessary function (ReadSOSCGTable) for input of 
previously-calculated S0(5) D S0(3) Clebsch-Gordan co- 
efficients. The function TruncateBasisMatrices may 
then be used to truncate the stored matrices of S0(3)- 
reduced matrix elements to lower t;i„ax, in order to re- 
duce the product space dimensions for ACM calcula- 
tions. The GcunmaHarmonic package also defines several 
functions (S05LSpaceSize, SOSBranching, SOSLSpace- 
Seniorities, etc.) based on the 80(5) D S0(3) dimen- 
sion and branching formulas of Appendix A (see the in- 
ternal program documentation). These functions are of 
use in indexing the SU(1, 1) x S0(5) basis functions and 
constructing the seniority contribution to the kinetic en- 
ergy operator. 



14 



5. CONCLUSION 

The computational methods described in this article, 

as implemented in the accompanying computer code, 
make possible the calculation of a sufficient set of S0(5) D 
S0(3) Clebsch-Gordan coefficients to support fully con- 
verged nuclear structure calculations with the algebraic 
collective model (ACM). The SU(1, 1) x S0(5) algebraic 
structure of the ACM basis permits matrix elements of 
an essentially unlimited set of potential and kinetic en- 
ergy operators to easily be constructed. The Bohr collec- 
tive model in the resulting calculational scheme is thus 
genuinely an algebraic collective model. The availability 
of S0(5) D S0(3) Clebsch-Gordan coefficients, in con- 
junction with a large body of analytic expressions for /3 
matrix elements [3, 4], makes it possible to algebraically 
construct the matrix elements for any collective model 
Hamiltonian expressible as a polynomial in the collective 
quadrupole moments and canonical momenta. Algebraic 
expressions for matrix elements of a variety of other oper- 
ators, including the term occurring in the Davidson 
potential [40, 41], have also been derived [3, 4]. 

The calculated S0(5) D S0(3) Clebsch-Gordan coeffi- 
cients now permit the diagonalization of the Bohr Hamil- 
tonian for potentials of essentially arbitrary 7 stiffness. 
This allows application of the ACM to the full range of 
nuclear quadrupole rotational- vibrational structure, from 
spherical oscillator to axial rotor to triaxial rotor. With 
an appropriate optimized choice of (3 basis functions [3], 
fully converged calculations can be carried out very effi- 
ciently, providing a valuable tool for studying collective 
motion in nuclei. 

The SU(1, 1) X S0(5) framework also opens the door for 
more detailed exploration of the formal relationship be- 
tween the Bohr collective model and the interacting bo- 
son model (IBM) [11], through the U(6) D U(5) D S0(5) 
and U(6) D S0(6) D S0(5) bases. Algebraic collective 
model methods have already been applied [13] to calcu- 
lation of collective states within the S0(6) limit of the 



IBM. The availability S0(5) D S0(3) Clebsch-Gordan 
coefficients enables such studies to be enhanced and ex- 
tended. 
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APPENDIX A: BRANCHING AND 
MULTIPLICITY FOR THE REDUCTION 
SO(5) D SO(3) 

The multiplicity of each S0(3) irrep (L) within the 
S0(5) irrep {v, 0) is needed in the code, both to determine 
the S0(5) D S0(3) branching for the basis and to convert 
between the {vaL) and {Li) labeling schemes. The result 
of Refs. [7, 42] can be expressed more compactly as 

dvL = (Li(^;-6)J +l)0.-b- ll{v-L + 2)\9,^L+2, (Al) 

where b = L/2 for L even or (L-|-3)/2 for L odd. The step 
function Ok is unity for > and zero otherwise. The 

dimension Dy^^^.^^i = Ylv^'o '^vL of the basis for a given 
L-space, truncated at seniority v<Vniax, is consequently 

-D^maxi = [/(^^max -b,3) + (t;i„ax - b + 1)] 

-/(w-i + 2,3)^,__i+2, (A2) 
where f{n,m) = J2k=ol^/''^\ given by 

/(n, m) = [n/mj [n+1- im( [n/mj + 1)]. (A3) 
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